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Abstract. A newly developed weak Galerkin method is proposed to solve parabolic equations. 
This method allows the usage of totally discontinuous functions in approximation space and pre- 
serves the energy conservation law. Both continuous and discontinuous time weak Galerkin finite 
element schemes are developed and analyzed. Optimal order error estimates in both and 
norms are established. Numerical tests are performed and reported. 



1. Introduction 
In this paper, we consider the initial-boundary value problem 

ut — • {aS/u) = f in ri, t & J, 
(1) u = g on d^l, t G J, 

n(-, 0) = ip in Q, 

where is a polygonal or polyhedral domain in {d = 2, 3) with Lipschitz-continuous boundary 
dQ, J = (0, i], and a = {aij{x,t))j_xd ^ [L°°{Q x J)]'^^ is a symmetric matrix-valued function 
satisfying the following property: there exists a constant a > such that 

For simplicity, we shall concentrate on two-dimensional problems only (i.e., d = 2). 

For any nonnegative integer m, let H"^{Q) be the standard Sobolov space, which is the collection 
of all real-valued functions defined on with square integrable derivatives up to order m with 
semi-norm 



I Jn 



H-- 
and norm 

m 

\mrn,n^{Y.\^\ln)'^. 

3=0 

For simplicity, we use || - || for the L?' norm. 

The parabolic problem can be written in a variational form as follows 

(2) {ut,v) + (aVu,Vv) = {f,v) Vv e H^{n), t e J, 

n(-,0) = ^, 

where u is called a weak solution if u € L'^(0,t; H^{Q,)) and ut € L^(0, t; and if u = g on 

an. 
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Parabolic problems have been treated by various numerical methods. For finite element methods, 
we refer to [1] and [2]. Discontinuous Galerkin finite element methods were studied in [3] and 0]. 
In [5] and [6], finite volume methods were presented, which were based on the integral conservation 
law rather than the differential equation. The integral conservation law was then enforced for small 
control volumes defined by the computational mesh. 

The goal of this paper is to consider a newly developed weak Galerkin (WG) finite element 
method for parabolic equation which is based on the definition of a discrete weak gradient operator 
proposed in [3. In this numerical method, the analysis can be done by using the framework of 
Galerkin methods, and totally discontinuous functions are allowed to be used as the approximation 
space. Furthermore, the approximation results also satisfy the energy conservation law. 

The rest of this paper is organized as follows. In Section 2, we introduce some notation and 
establish a continuous time and discontinuous time weak Galerkin finite element scheme for the 
parabolic initial boundary- value problem ([T]). In Section 3, we prove the energy conservation law 
of the weak Galerkin approximation. Optimal error estimates in both and norms are proved 
in Section 4. The paper is concluded with some numerical experiments to illustrate the theoretical 
analysis in Section 5. 

2. The Weak Galerkin Method 

In this section we design a continuous time and a discontinuous time weak Galerkin finite element 
scheme for the initial-boundary value problem ([T]). We consider the space of discrete weak functions 
and the discrete weak operator introduced in [7]. Let 7/i be a quasiuniform family of triangulations 
of a plane domain Q and T be each triangle element with as its interior and dT as its boundary. 
For each T G Th, let Pj{T^) be the set of polynomials on with degree no more than j, and 
Pi{dT) be the set of polynomials on dT with degree no more than /, respectively. Let PjiT) be 
the set of homogeneous polynomials on T with degree j. The weak finite element space Sh{j,l) is 
defined by 

Sh{j,l) := {v = {vo,Vb} : vo € P,(r°), Vb G Pi{e) for all edge e C 9T, T G %}■ 

Denote by S^{j, I) the subspace of Sh{j, I) with vanishing boundary value on dQ; i.e., 

ShU^l) '■= {v = {vo,Vb} G Sh{j,l), Vblarnan = for all T G Th}. 

Let Eh = {q G (L'^i^))'^ ■ q|T e V{T,r) for all T G %}, where V{T,r) is a subspace of the set 
of vector-valued polynomials of degree no more than r on T. For each v = {v(),Vb} G Sh{j,l), the 
discrete weak gradient V^^f G E/i ^ each element T is given by the following equation: 

(3) / VdV ■ qdT = - VoV-qdT+ / v^ci-nds, VqGy(r,r). 
Jt Jt JdT 

It is easy to see that this discrete weak gradient is well-defined. 

To investigate the approximation properties of the discrete weak spaces Sh{j, I) and 
three projections in this paper. The first two are local projections defined on each triangle T: one 
is QhU = {Qqu, Qbu}, the projection of H^{T) onto Pj(T^) x Pi{dT) and another is Rh, the 
projection of [L^(T)]^ onto V{T,r). The third projection li^ is assumed to exist and satisfy the 
following property: For q G H{div, ft) with mildly added regularity, Il/iq G H{div, Q) such that 
Il/iq G V{T, r) on each T G Th, and 

(V • ci,vo)t = (V • nhq,vo)T, yvo G Pj(rO). 

It is easy to see the following two useful identities: 

(4) Vd{Qhu) = Rh{Vu), yueH\T), 
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and for any q G H{div, 0) 

(5) ■ ^o)^ = E (n^q, Vdr;)T, Vr; = {t-o, Vb} G Z). 

The discrete weak spaces Sh{j,l) and need to possess some good approximation properties 
in order to provide an acceptable finite element scheme. In [7j, the following two criteria were given 
as a general guideline for their construction: 

PI: For any v € Sh{j,l), if ^d'f = on T, then one must have v = constant on T; i.e., 

^^0 = ^^fe = constant on T. 
P2: For any u € Hq(Q) n H'^'^^{Q), where < m < j + 1, the discrete weak gradient of the 
projection Q^u of u in the discrete weak space Shij,l) provides a good approximation 

of Vu; i.e., \\S/(iiQhu) — Vu|| < C/i™ ||u||m+i holds true. 

Examples of Sh{j,i) and satisfying the conditions PI and P2 can be found in [7]. For 
example, with V{T,r = j + 1) = [Pj{T)]'^ + Pj(T)x being taken as the usual Raviart-Thomas 
element [U] of order j, one may take equal-order elements of order j for Sh{j, the interior and the 
boundary of each element T. The key of using the Raviart-Thomas element for V{T, r) is to ensure 
the condition PI. The condition P2 was satisfied by an operator equation RhS = ^dQh which 
has been established in [7]. It was shown later in [H] that the condition PI can be circumvented 
by a suitable stabilization term. Consequently, the selection of V{T^r) and Sh{j,l) becomes very 
flexible and robust in practical computation. It even allows the use of finite elements of arbitrary 
shape. 

The main idea of the weak Galerkin method is to use the discrete weak space Sh{j,l) as testing 
and trial spaces and replace the classical gradient operator by the weak gradient operator in 

©• 

First, we pose the continuous time weak Galerkin finite element method, based on the weak 
Galerkin operator ([3]) and weak formulation ([2]), which is to find Uh{t) = {uo{-,t),Ub{-,t)}, belonging 
to Sh{j,l) for t > 0, satisfying Ub = Qtg on d^l,t > and Uh{0) = Qhi^ iii and the following 
equation 

(6) {{uh)t,vo) +a{uh,v) = {f,vo) = {wq, S /), t > 0, 
where a(-, •) is the weak bilinear form defined by 

a{w,v) = {aVdW,S/dv), 
which is assumed to be bounded and coercive, i.e., for constant a,/3,7 > 

\a{u,v)\<p\\Vdu\\\\Vdvl 
a{u,u) > Q||Vdu||^, 

and that 

\{atVdU,Vdv)\ <7l|Vd^^||||Vrfi;||. 

We now turn our attention to some discrete time Weak Galerkin procedures. We introduce 
a time step k and the time levels t = tn = nk, where n is a nonnegative integer, and denote 
by [/"■ = UJ^ S Sh{jJ) the approximation of u{tn) to be determined. The backward Euler Weak 
Galerkin method is defined by replacing the time derivative in ([6]) by a backward difference quotient, 
or, if BU"^ = - U''-^)/k, 

(7) (5C/",^;o) + a(C/",^;) = {f{tn),vo) Vi; = {vo,Vb} G Sh{j,l), n > 1, ^7° = Q^V, 
i.e., 

{U^,vo) + kaiU'', v) = ([/"-I + kfitn),vo), yv = {vo,Vb} G ShU, I). 
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3. Energy Conservation of Weak Galerkin 



This section investigates the energy conservation property of the weak Galerkin finite element 
approximation u/j. The increase in internal energy in a small spatial region of the material i.e., 
control volume, over the time period \t — At, t + At] is given by 

r r rt+At i- 

I u{x,t -\- At) dx — / u{x , t — At) dx = / / utdxdt. 

Jk Jk Jt-At Jk 

Suppose that a body obeys the heat equation and, in addition, generates its own heat per unit 
volume at a rate given by a known function / varying in space and time, the change in internal 
energy in K is accounted for by the flux of heat across the boundaries together with the source 
energy. By Fourier's law, this is 

t+At j- rt+At j- 

I q-ndsdt+ fdxdt, 
t-At JdK Jt~At JK 

where q = —aS/u is the flow rate of heat energy. The solution u of the problem ([1]) yields the 
following integral form of energy conservation: 

i-t+At r- rt+At I- rt+At i- 

(8) / j Utdxdt + j q ■ ndsdt = f dxdt 



it-At Jk Jt-At JdK Jt-At Jk 

where the Green's formula was used. We claim that the numerical approximation from the weak 
Galerkin finitel element method for ([1]) retains the energy conservation property 

In ([6]), we chose a test function v = {vo,Vb = 0} so that = 1 on i(' and vq = elsewhere. After 
integrating over the time period , we have 

rt+At j- rt+At r rt+At r 

(9) / I Utdxdt + j I aVduVdV dxdt = f dx dt. 

Jt-At Jk Jt-At Jk Jt-At Jk 

Using the definition of operator Rh and of the weak gradient in ([3|), we arrive at 

/ aVduVdvdx= / Rh{aVdUh) -Vdvdx = - V ■ Rh[aV dUh) dx = - I Rh{aV dUh) ■ n ds . 
Jk Jk Jk JdK 

Then by substituting in (l9|), we have the energy conservation property 

rt+At j- rt+At r rt+At r 

j I Utdxdt + j / —Rh{aS/ dUh) ■ ndsdt = j / fdxdt, 
Jt-At Jk Jt-At JdK Jt-At Jk 

which provides a numerical flux given by 

qh-n = -RhiaVdUh) ■ n. 

The numerical flux qh ■ n is continuous across the edge of each element T, which can be verified by 
a selection of the test function v = {vQ,Vf,} so that f o = and Vf, arbitrary. 



4. Error Analysis 

In this section, we derive some error estimates for both continuous and discrete time weak 
Galerkin finite element methods. The difference between the weak Galerkin finite element approx- 
imation Ufi and the projection QhU of the exact solution u is measured. We first state a result 
concerning an approximation property as follows. 

Lamma 4.1. For u € H^^^{Q) with r > 0, we have 

\\Uh{aVu) -aRhiVu)\\ < C h'' \\u\\i+r ■ 
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Proof. Since from (jH) we have i?/i(Vti) = Vd{Qhu), then 

\\Uh{aVu) - aRhiVu)\\ = ||n,,(aVu) - aVd{Qhu)\\. 

Using the triangle inequahty, the definition of Hh and the second condition P2 on the discrete weak 
space Sh{j,l), we have 

||n/j(aVn) — aV(i{Qhu)\\ < ||n/i(aVu) — aVn|| + \\aVu — a'S/d{Qhu)\\ < C/i^||n||i+j.. 

□ 

4.1. Continuous Time Weak Galerkin Finite Element Method. Our aim is to prove the 
following estimate for the error for the semidiscrete solution. 

Theorem 4.1. Let u € H'^^'''{yL) and be the solutions of (JD and respectively. Denote by 
e := Uh — QhU the difference between the weak Galerkin approximation and the projection of the 
exact solution u. Then there exists a constant C such that 

||e(-,t)f + [ a\\Vdefds<\\e{-,0)f + Ch^'' [ \\u\\l+^ds, 
Jo Jo 

^Vtf ds + fl|V,e(-,i)f + (1 + 



and 



< a||Vrfe(-,0)f + (l + ^)||e(-,0)f +C/i2-(|ln(.,0)||?+, + ||n||2^, + ^ ||n||?+,ds + ^ htfi+rds). 

Proof. Let v = {vo,Vb} G S'^(j, /) be the testing function. By testing (HD against vq, together with 
RhCVu) = VdiQhu) for u G H^, where R is the local projection onto V{T,r), and {Qout^vo) = 
{ut, Vq), we obtain 

(/,^^o) = {ut,vo) + "^{-V ■ {aVu),vo)T 
= {ut,vo) + (Uh{aVu),S/dv) 

= (QhUt, vo) + {Uh{aVu) - aRh{Vu),Vdv) + {aV d{Qhu),V dv) 
Since the numerical solution also satisfies the heat equation, we have 

(/, ^^o) = {{Uh)t, Vq) + a{uh,v). 

Combining the above two equations we obtain 

(10) {{uh - Qhu)t, Vq) + a{uh - QhU, v) = (n/j(aVn) - aRh(yu),S/dv), 

which shall be called the error equation for the WG for the heat equation. 
Let e = Uh — QhU be the error. Use t; = e in the error equation, we obtain 

(et, e) + a(e, e) = (Uh{aVu) - aRh{Vu),S/de)- 

By Cauchy-Schwarz inequality and the coercivity of the bilinear form, we have 



It follows that 



Ij^W^f + «l|Vdef < ^\\Uh{aVu) - aRh{Vu)f + |||Vrfef . 

^||e|p+a||Vdef < -\\UhiaVu) - aRh(yu)f , 
at a 



and hence, after integration, 

(11) \\ef+ ['a\\Vdefds<\\e{;0)f + - [' \\nh{aVu) - aRh{Vu)f dt. 

Jo a Jo 
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Then by Lemma l4.ll we have the assertion. 

In order to estimate V^e, we use the error equation with v = {uh — Qhu)t = st- We obtain 

{et, et) + a(e, et) = {Ilh{aVu) - aRhiVu), VdCt) 
= ^(UhiaS/u) - aRh{Vu),Vde) - (UhiaVut) - aRh{Vut), V^e) - (Uh{atVu) - atRh{Vu),Vde). 

By Cauchy-Schwarz inequahty, this shows that 

\\etf + ^i-^a{e,e) - {atVde,Vde)) < ^{Uh{aVu) - aRhiVu),Vde) 

+ ^\\Uh{aVut) - aRh{Vut)f + ^HV^ef 
2a 2 

+ J-lln^(atVn) - atRh{Vu)f + ^HV^ef , 
2a 2 

i.e., 

1 1 d 

||etf + -^a(e,e) < - (at V^e, V^e) + —(11,, (aVu) - ai?/,(Vu), V^e) 

+ ^\\Uh{aVut) - ai2ft(Vnt)f + ^llV^ef 
2a 2 

+ ^WUhiatVu) - atRh{'^u)f + ^\\\7def. 
2a 2 

Thus, integrating with respect to t and together with the coercivity and boundedness yields 



f\\etfds + ^\\Vde{-M 



< fl|Vde(,0)f 



+ {UhiaVu{;t) - aRhiVu{;t)), Vde{;t))-{UhiaVu{; 0) - aRh{Vu{;0)),Vde{;0)) 

+ 7^1 WllhiaVut) - aRh{Vut)f ds + ^ [ \\Uhiat\/u) - atRh{'^u)f ds 
2a Jq 2a Jq 

+ + W'^defds, 
By adding {a + j/2)/a = 1 + 2^ times inequahty (fTT]l to the above inequahty, we arrive at 



f\\et\\'ds + ^\\Vdei;t)f + ii + 
Jo ^ 

< fllVde(-,0)f + (l + ^)||e(.,0)f 

+ {Uh{aVu{; t)) - aRh{Vu{-,t)), Vde{-,t))-{Uh{aVu{-, 0)) - aRh{Vu{-, 0)), Vde(-, 0)) 

+ 77- / \\Ilh{aVut) - aRh{Vut)f ds + ^ [ \\Uh{atVu) - atRh{Vu)\\^ ds 
2a Jo 2a Jq 

+ (- + Tr^)/ ¥^h{aSIu)-aRh{Vu)fdt, 
a 2a^ Jq 



Next, we use Cauchy-Schwarz inequality to obtain 

l\\e,fds + ^\\V,e{;t)f + {^ + ^)\\ef 
< /3||V,e(,0)f + (l + ^)||e(,0)f 

+ -\\UhiaVu{;t)) - aRhiVui-,t))f + ^||n,,(aV^x(-, 0)) - ai?,,(Vu(-, 0)) f 

+ ^ [ \\Uh{aVut) - aRh{Vut)f ds + ^ [ \\Uh{atVu) - atRh{Vu)f ds 
2a Jo 2a Jq 

+ (- + A) / \\Uh{a\7u) - aRh{Vu)f dt. 
a 2a^ Jq 

Then by Lemma l4.ll the proof is complete. □ 

4.2. Discrete Time Weak Galerkin Finite Element Method. We begin with the following 
lemma which provides the Poincare inequality with the weak gradient operator. 

Lamma 4.2. Assume that (p G S^{j,j), then there exists a constant C such that 

|2 / r^llw _J,ll2 



122 < c||v, 

where Vacj) G ViT, r = j + 1) = [Pj(r)]2 + P, (T)x. 

Proof. Let cp^ G S^{j,j), a piecewise constant, be the cell average of (p. Let cpi G S'jl{j,j) dHQ^Q) 
be a continuous piecewise polynomial with vanishing boundary trace lifted from a discontinuous 
piecewise polynomial as follows. Let Gj{T) be the set of all Lagrangian nodal points for Pj{T). At 
all internal Lagrangian nodal points xj G Gj{T), we set (pi{xj) = (fP{xj). At boundary Lagrangian 
points Xj G Gj{T) PI dT, we let (f)i[xj) be either a trace of (jP from any side. At global Lagrangian 
points Xj G dT n 50, we set (j)i{xj) = 0. Let denote the jump on the edge e, i.e., [0]^ = 
(p^lTini + (/>'^|t2?^2 = {(P^Iti — (p^)ni + {(I>^\t2 ~ (l>^)'n'2- By the classical poincare inequality for cpi, we 
have 

I2 < U-Mh + Ui\\h 



From Lemma 4.3 in [10], we have 



T'€TiT) eee{T) 

where T{T) denotes the set of all triangles in T having a nonempty intersect with T, including 
T itself and £{T) denotes the set of all edges having a nonempty intersection with T. Note the 
elementary fact that 



T 



where Xi and Xj are two end points of an edge. Combine above inequalities, we obtain 

(12) Wl^ <f^(E^^ii^'^°iiT+E^^iiMiie+E^^"'r'""' 



T e e 

To this end, we want to bound all three terms on the right hand side by ||Vd 
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Let us recall that the definition of the weak gradient V d(t> 

/ Vd(p-qdx = - (j)°V-qdx+ / • nds, Vq G T/(r, j + 1), 

Jt Jt Jot 

and by using integration by parts, we have 

(13) [ \/d<l)-(ldx= [ -qdx- [ ((/.° - /)q • Jids, VqGy(r,j + l). 

Jt Jt JdT 

In order to bound ||V(/)^|| and ||[0]||e by HVdi^ll, we let q in ()13p satisfy the following 

/ q-rdx = [ Wfp^-rdx Vr € [Pj(r)]2, 
Jt Jt 

[ q-nfids = -h~^ f {(jP - (l)^)^^Lds G Pj+i[9T], 
JdT JdT 

where by Lemma 4.1 in [11], q and (f) have the norm equivalence of 



(14) l|q||L2{T) ~ l|v/||L2(r) +/i ^11/ -/||^2(ar)- 

Then from (jlSp . we have 



f Vd(t>-qdx = [ V4>°-qdx- [ {(1)^ - (t>^)q ■ nds 
Jt Jt JdT 



■ V/ dx + / (/ - /)2 ds 

T JdT 

= \\v^YT + h-V-A\dT- 

Using (fTH]) and i^^., we have 

||Vd<A||T>^(V,</.,q)T = ^(||V</.°||^ + /i-i/-0lir) 

> c{\\vcpYT + h-V-A\oT), 

then together with (I12p . we have the assertion. □ 

With the results established in Lemma 14. 2^ we are ready to derive an error estimate for the 
discrete time weak Galerkin approximation Uh as the following theorem. 

Theorem 4.2. Let u G //^"^""(il) and f/" be the solutions of ([7]j and respectively. Denote by 
:= [/" — Qhu{tn) the difference between the backward Euler weak Galerkin approximation and 
the projection of the exact solution u. Then there exists a constant C such that 

||e"f + VaA;||Vde^'f < |le°f + C(/i2'^||m||?+^ + A;2 / \\uuf ds), forn>0, 

and ^ 

l|Vde"f < C (lleOf + llVrfgOf + h^'{\\u\\l+, + \\ut\\l+,) + Jj hufds^ , 

where \\u\\i^r = max {||ii(ij)||i+r} ctnd \\ut\\i+r = max {||uj(tj)||i+r}. 



Proof. A calculation corresponding to (|TO]l yields 
(15) 

{d{U''-Qhu{tn)),vo)+a{U''-Qhu{tn)),v) = {ut{tn)-du{tn),vo)+{Uh{aVu{tn))-aRh{Vu{tn)),Vdv), 
i.e., 

(16) (9e", vo) + a(e", t;) = (nt(t„) - du{tn),vo) + {Uh{aVu{tn)) - aRh{Vu{tn)),Vdv). 



Let Wi = ut{tn) — du{tn), and W2 = n/i(aVu(t„)) — ai?/i(Vn(t„)), and choosing v = e"" in ([16 
we have 

(ae^ e") + a(e", e") = K, e") + K, V,e"). 
By coercivity of the bihnear form and Cauchy-Schwarz inequahty, we obtain 

||e"f - (e"-\e") + ak\\Vy'f < A;|K||||e"|| + k\\w'^\\\\Vde^, 



I.e. 



e"f + aA;||V,e"f < ^lle"-^^ ^ Kn^^2 ^ A:||^«?||||e"|| + A;||u;^||||V,e"||. 



2 2 
By the Poincare inequahty of Lemma 14.21 we have 



Then by triangle inequahty, it follows 



so that 



IWelf + ^l|V.e"f < l\\e-~r + ^(^^H^i II + 11^2 11)^- 



a a 



and, by repeated application, 

n 

(17) ||e"f + ^afcllV.e^-f < He^f + — ^ H^l f + " E 



i\ X ^ II 1 1 1 o > II J 1 1 2 

i/^Q I 

a ■' — ' ■ - ■ Q, — i ■■ 



We write 

(18) km[ = kut{tj) - {u{tj) - u{tj_i)) = / (s - ij_i)'Ujt(s) ds. 



I.e., 



i-i 



SO that 



I'^ill^ = i ^\ {s -'tj~i)utt{s)dsY dx 



(19) < j j {s-tj-ifds I uit{s)dsdx 



Substitute p9p into p7p and together with Lemma |4.H we have the error estimate for ||e"||. 

In order to show an estimate for Ve" we may choose instead v = Be"" in error equation ()17p to 
obtain the following error equation 

(ae*^, 5e") + o(e", ae") = {w^^, Oe") + {w^, Vd^e"). 

The second term on the right hand side can be written as 



Then the error equation becomes 

fepe^f + (aVrfe'^,Vrfe") = (aV^e", V^e"-!) + fe^, 5e") 

By triangle inequahty, we have 

+ + ^l|V.e"-if, 

and, after cancelation and by repeated application, 
i(aV,e",V,e") < i(aVrfe°, V^e") 

+ f Ell^ill' + (^2,V,e")-K0,Vrfe°) 

j=l j=l j=l 



which is 



(20) + \ jZ\Hf + + f l|V.e"f + hwlf + f ||V.e°f 

+ f E 11(^2)* - dwif + 1 E + ^E 

j=l j=l j=l 

followed by triangle inequality and boundness and coercivity of the bilinear form. By the similar 
process as in (|19p and Lemma l4.ll we have 

- Bwif < Ck r \\{w2)ufds < Ckh^' r \\utt\\l+rds. 

Then by substituting (|19p . (jl7p and the above inequality into (|20p . we have the error estimate 
for llVrfe"!! as the following 

l|Vde"f < C [\\eY + llV^eOf + + \\ut\\l+r) + k"" j'j \\uufds^ , 

which completes the proof. □ 

4.3. Optimal Order of Error Estimation in L?. The optimal order of error estimation for 
||Vde|| and llV^e"!! was obtained in section 4.2. In order to get an optimal order of error estimate 
in L^, the idea, similar to Wheeler's projection as in |12] and [1], is used where we define an elliptic 
projection Eh onto the discrete weak space Sh{j,l) as the following: 

(21) {aS/dEhV,S/dX) = {-'^-{a\/v),x), ^X^ShUJ), fove^veH^. 
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In view of the weak formulation of the elliptic problem, 

(22) - V • (aVv) = F in Q, 

u = g on 50, 

this definition may be expressed by saying that EhV is the weak Galerkin finite element approx- 
imation of the solution of the corresponding elliptic problem with exact solution v. By the error 
estimate results in [7], we have the following error estimate for EhV. 

Lamma 4.3. Assume that problem has the H^^^ regularity (s G (0,1])- Let v G H^^^' he 
the exact solution of 112^) . and E^v he a weak Galerkin approximation of v defined in i21\) . Let 
QhV = {Qov,Qbv} be the projection of v in the corresponding finite element space. Then there 
exists a constant C such that 

(23) IIQo^^ - Ehv\\ < C{h'+'\\F - QoF\\ + h^+^\v\\r+i), 
and 

(24) W'^diQhV-EhV^ <Ch'-\\v\\r+i. 

Throughout this section the error in the parabolic problem ([T]) is written as a sum of two terms, 

(25) Uh{t) - Qhu{t) = 9{t) + p{t), where 9 = Uh - EhU, p = EhU-QhU, 

which will be bounded separately. Notice that the second term is the error in an elliptic problem 
and then can be handled by applying the results in Lemma 14.31 Then our main goal here is to 
bound the first term 6. 

Following the above strategy, the error estimate for continuous time weak Galerkin finite element 
method in and are provided in the next two theorems. 

Theorem 4.3. Under the assumption of Theorem \4.1\ and the assumption that the corresponding 
elliptic prohlem has the H^^^ regularity (s S (0,1]), there exists a constant C such that 

\\uh{t) - Qhu{t)\\ < \\uh{^)-Qhu{m+Ch'^'{\Mr+i+ f \\ut\\r+ids) 

Jo 

+ C/i^+i(||/(0) - Qo/(0)|| + \\ut{0) - Qoutmi) 
+ Ch'^'[l {\\ft-Qoft\\ + \\uu-QoUtt\\)dsY 

Proof. We write the error according to (I25p and obtain the error bound for p easily by Lemma 14.31 
as the following 

(26) IIpII < C{h'+\\\f - Qo/ll + lint - QoUtW) + /i'-+^(||^||.+i + /* |lnt||,+i ds)). 

Jo 

In order to estimate 9, we note that by our definitions 

{Ot,x) + a{9,x) = {uh,t,x) + a{uh,x) - {EhUt,x) - a{EhU,x) 

= (.f,x) - {EhUt,x) - a{EhU,x) 

= if,x) + iy ■iaVu),x)-iEhUt,x) 

= {ut,x) - {EhUt,x) 

= {QhUt,x) - {EhUt,x) 

= -{pt,x), 

which is 

(27) {0t,x) + a{9,x) = -{Pt,x), Vx S 5,,(i, /), t > 0, 
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where we have used the fact that the operator Eh commutes with time differentiation. Since 
9 S Sh{j, I), we may choose x = ^ hi (|27p and obtain 

i9t,e) + ai9,9) = -{pt,9), t>0. 

Since 

a{0,0) > a\\V9f > 0, 

we have 

i^||^||2 = ||0||^||^||<|| llll^ll 

and hence ^ 

||e(i)ll< 11^(0)11+ [ \\pt\\ds. 

Jo 

Using Lemma 14.31 we find 

\\9{0)\\ = \\uhiO) - Ehu{0)\\ = \\uh{0) - Qhumi + \\Ehu{0) - Q/,n(0)|| 

(28) < ||n^(0) - Q^^z(0)|| + C[h'+\\\f{0) - Qo/(0)|| + ||nt(0) - Qont(0)||) + 
and since 

(29) \\pt\\ = \\EhUt - QhUtW < C[h'+\\\ft - QoftW + \\uu - QoUu\\) + h'^+'WutWr+i], 

the desired bound for ||0(t)|| now follows. □ 

Theorem 4.4. Under the assumption of Theorem \4.3\ and the assumption that the coefficient matrix 
a in is independent of time t, there exists a constant C such that 

W^diuhit) - Qhu{t))f < 4/3||Vrf(tx^(0) - Qhumf + C/i^^XllV'll'+i + \\ufr+i) 

+ C7/i2(^+i) / (\\f^ - Qof^W + \\uu - QoUttW? ds + Ch^^^'+'^ [ \\uS_^.^ds. 
Jo Jo 

Proof. As in the proof of Theorem 14.31 we write the error in the form (j25p . Here by Lemma 14.31 

(30) \\VdPm<Ch'\\u\\r+i. 
In order to estimate V^^, we may choose x = ^t- We obtain 

{9t,9t) + a{9,9t) = -{pt,9t). 
Since coefficient matrix a in the bilinear form a(-, •) is independent of time t, we have 



Id, 1„ „o 1, 

I /J 



11^*11' + 7T3T(«Vde,Vrf0) = -{pt,9t) < -\\ptf + -1'^ 



2dt" " ' - ' u - 2 
so that 

j^{aVd9,Vd9) < \\ptf. 

Then by integrating with respect to time t and using the coercivity and boundedness of the bilinear 
form, we obtain 

a\\Vd9f < {aVde,Vde)<{aVd9{0),Vd9m+ f \\ptf ds < (3\\Vde{0)f + f \\pt\? ds 

Jo Jo 

< m\'^d{uh{0)-Qhu{0))\\ + \\Vd{Ehu{0)-Qhu{0))\\f + [ \\ptf ds. 

Jo 

Hence, in view of Lemma 14.31 and (j29p . we have 

lr+1 



\Wd9it)f < 2/3||Vrf(n;,(0)-Q,n(0))f + C/i2-n„,,ii2 



+Ch^{s+i) f\y^_Q^f^\\ + \\u^^_Q^u^^\\fds + Ch'<^^+''> [ 
Jo Jo 



t 

utWlj^i ds, 
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which completes the proof. □ 

Next, we prove the error estimate for the backward Euler weak Galerkin method. 

Theorem 4.5. Let u S H^^^'{Q) and [/" be the solutions of ([7]j and respectively. And let QhU 
be the projection of the exact solution u. Then there exists a constant C such that 

W-Qhu{tn)\\ < \\U^ -Qhum\+Ch^^"{U\\r+i+ r \Wt\Uids) 



< 


||C/°- 


+ 


Ch'+\ 


+ 


Ch'+\ 


+ 


Ch'+^ 






+ 


Ck / 




Jo 







\uu\\ ds. 

Jo 

Proof. In analogy with (j25p . we write 

f/" - Qhu{tn) = - Ehu{tn)) + {Ehu{tn) - Qhu{tn)) = r + p", 
where p" = p{tn) is bounded as shown in (j26p . In order to bound 0", we use 

(5r,x) + a(r,x) = {dU^,x)+a{U^,x)-{dEHu{tn),x)-a{E,,u{tn),x) 
= {f{tn),x) - {dEhu{tn), x) - a{Ehu{tn), x) 
= (/(tn), X) + (V • {aVu{tn)),x) - {dEhu{tn),x) 

= {ut{tn),X) - {dEhU{tn), X) 

= {Ut{tn) - du{tn),x) + {du{tn) - dEhu{tn),x), 

i.e., 

(31) (5r,x) + a(r,x) = K,x), 

where 

W'' = (Utitn) - du{tn)) + {du{tn) " dEhU{tn)) = w'^ + . 

By choosing x = in (|3ip and the coercivity of the bilinear form, we have 
or 

||^n||2 _ (gn-l^^n) < A; || U;" || || 0" || , 

SO that 

||^"|| < +A;||u;"|[, 
and, by repeated application, it follows 

n n n 

ll^ll < |l^°|| + A:^||u;^'|| < \\0^\\ + kJ2\\w{\\ + kJ2\\wi\\- 
j=i j=i j=i 

As in (j28p . 6^ = 6{0) is bounded as desired. By using the representation in (|18p . we obtain 

" rtj rtn 

^Xl ll^ill - X] II / {s -tj_i)utt{s)ds\\ < k \\utt\\ds. 
-•1 Jt,_i Jo 



j=i j=i •'^i-i 



We write 



(32) wi = du{tn) - Ehdu{tn) = {Eh - I)k~^ / ' utds = k-^ T [En - I)ut ds, 
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and, by ([29]) we have 

n n „4 . 

fc^ll^ll < I C[h'^\\\ft-Qoft\\ + \Wtt-QoUtt\\) + h'+'\\ut\\r+l]ds 

i=i j=i 

fin f'tn 

= C[h'+' (.\\ft-Qoft\\ + \\uu-QoUu\\)ds + h^'+' \\ut\\r+ids], 
Jo Jo 

Thus, together with the estimate of p in ()26p . we have the assertion. □ 

Theorem 4.6. Under the assumption of Theorem \4-5\ and the assumption that the coefficient 
matrix a is independent of time t, there exists a constant C such that 

||V,([/" - Qhu{tn))\? < 2||V,(C/0 - Qhumf + Ch^'^iMl^, + ||n||,V) 

+ / (\\f^ - Q,f^\\ + \\uu - QoUttWf ds + Ch^'-^+'^ \\utfr+ids 



10 JO 

j-tn 

Ck'^ / \\uu\\'^ ds. 
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Proof. In order to estimate V^^", we choose x = iii (|3ip . and it is easily seen that 

aailVrf^f < \\w 



so that 

an II 2 ^ II „,.n ||2 



By repeating the apphcation, we have 

l|Vrf0"f < wv^eY + -y Ik^'f < liv<i^°f + 2- y + 2-y 



j=0 j=0 3=0 

As in (|32]). we obtain 



k\\w{f = k I {k~^ j ptdsfdx< I [j pfds)dx< j \\pt\? ds. 
Together with ()19p . (j29p and (jSOp . we have the assertion. □ 

5. Numerical Experiments 

In section 2, we mentioned that the discrete weak space Shij,l) and in the weak Galerkin 
method need to satisfy two conditions. In [7], they propose several possible combinations of Sh{j, I) 
and ^^j. Through out this section we use a uniform triangular mesh 7h, the discrete weak space 
Sh{0, 0), i.e., space consisting of piecewise constants on the triangles and edges respectively, and 
with V{T, 1) to be the lowest order Raviart-Thomas element RTq{T) in the weak Galerkin method, 
which were used in [13] for the numerical studies of the weak Galerkin method for second-order 
elliptic problems. We also adopt the various norms used in |13] to present the numerical results of 
the error e^ between the projection of the exact solution QhU and the numerical solution u/,. 

Example 1. As the first example, we consider the following heat equation in = (0, 1) x (0, 1), 

ut-V ■ (aVu) = f in fi, for t > 0, 
(33) u = g on dft, for t > 0, 

u(-, 0) = ip in Q, 
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where a 



1 
1 



and /, g and -0 are determined by setting the exact solution u = sin(27r(t + 

1) + 7r/2) sin(27rx + 7r/2) sin(27ry + 7r/2). 

For this inhomogeneous Dirichlet boundary condition, with a uniform triangular mesh Th, we 
chose the approximation space 

r e Pom, for all T € %, 

Sh = lv = {vo, Vb} : Vb G Po(e) for all T e % and e C dT <^ dn, 
[ Vb = gh for all T G T; and dT G dQ 

where gh is the interpolation of g in the space Po{dQ,). In the test, k = h and k = h? are used to 
check the order of convergency with respect to time step size k and mesh size h respectively, since 
the convergence rate of the error is dominated by that of the time step size k when k = h, and the 
convergence rate of the error is dominated by that of the mesh size h when k = h? . The results 
are shown in Table 1 and Table 2. Since the exact solution is smooth, we observe the optimal 
convergence rates in both and weak Galerkin norms for the Dirichlet boundary data type 
initial boundary value problem, which is consistent with the theoretical results shown in section 4 
and 5. 

Table 1 . Convergence rate for heat equation with inhomogeneous Dirichlet bound- 
ary condition with k = h 



h 


\\^h\\[oo,T} 


\\^h\\{oo,dT} 


l|Vde,,|| 


N^/i {L2,T} 


{L2,9T} 


1/8 


1.38e-01 


1.44e-01 


2.41e-01 


4.90e-02 


8.78e-02 


1/16 


6.97e-02 


7.26e-02 


8.34e-02 


2.20e-02 


4.03e-02 


1/32 


3.47e-02 


3.56e-02 


2.97e-02 


1.05e-02 


1.94e-02 


1/64 


1.72e-02 


1.75e-02 


l.lOe-02 


5.16e-03 


9.54e-03 


1/128 


8.59e-03 


8.65e-03 


4.27e-03 


2.56e-03 


4.74e-03 


0{K')r = 


1.0012 


1.0138 


1.4550 


1.0643 


1.0533 



Table 2. Convergence rate for heat equation with inhomogeneous Dirichlet bound- 
ary condition with k = h'^ 



h 


h^?* {oo,T} 


k/i {oo.aT} 




N^/i {L2,T} 


^/i {L2,9T} 


1/8 


7.11e-02 


8.30e-02 


7.81e-02 


3.28e-02 


5.39e-02 


1/16 


1.89e-02 


2.28e-02 


1.84e-02 


8.53e-03 


1.38e-02 


1/32 


4.79e-03 


5.84e-03 


4.52e-03 


2.16e-03 


3.49e-03 


1/64 


1.21e-03 


1.48e-03 


1.13e-03 


5.43e-04 


8.79e-04 


1/128 


3.64e-04 


4.31e-04 


2.88e-04 


1.49e-04 


2.47e-04 


0{h'')r = 


1.9025 


1.8994 


2.0213 


1.9454 


1.9418 



Next, we consider this heat problem with a mixed boundary condition 

u = g on dQi, 
Vii ■ n + u = 0, on di}2, 

where dQi R = 05 and d^i U 80,2 = 90,. The exact solution is set to be n = sin(27r(t^ + 
1) + 7r/2) sin(7ry)e~^, where dil,2 is the boundary segment x = 1 and dOi is the union of all other 
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Table 3. Convergence rate for heat equation with Robin boundary condition with 
k = h 



h 


{oo,T} 


^/i {oo,dT} 






^/i {L2,9T} 


1/8 
1/16 
1/32 
1/64 
1/128 


1.65e-01 
l.OOe-01 
5.54e-02 
2.92e-02 
1.50e-02 


1.72e-01 
1.02e-01 
5.59e-02 
2.93e-02 
1.50e-02 


1.80e-01 
8.55e-02 
4.01e-02 
1.91e-02 
9.24e-03 


9.86e-02 
5.86e-02 
3.22e-02 
1.69e-02 
8.67e-03 


1.83e-01 
1.09e-01 
5.95e-02 
3.13e-02 
1.60e-02 


0{K')r = 


0.8656 


0.8793 


1.0713 


0.8767 


0.8789 



Table 4. Convergence rate for heat equation with Robin boundary condition with 



h 


'^h {oo,T} 


{oo.dT} 


^deh 




'^h {L2,9T} 


1/8 
1/16 
1/32 
1/64 
1/128 


3.18e-02 
8.29e-03 
2.10e-03 
5.24e-04 
1.39e-04 


3.90e-02 
l.Ole-03 
2.54e-03 
6.34e-04 
1.62e-04 


2.61e-02 
6.28e-03 
1.55e-03 
3.88e-04 
1.06e-04 


1.88e-02 
4.87e-03 
1.23e-03 
3.08e-04 
8.79e-05 


3.57e-02 
9.22e-03 
2.32e-03 
5.83e-04 
1.65e-04 


0{K')r = 


1.9591 


1.9785 


1.9875 


1.9352 


1.9383 



boundary segments. For the mixed boundary data type of initial boundary value problem, we also 
achieved the optimal convergence rates of the error in all norms as shown in Table 3 and Table 4. 

Example 2. For the second example, we consider the parabolic problem (j33p of full tensor with 



Dirichlet boundary condition, where coefficient matrix a 



+ + 1 



xy 



which 



xy X" + y" + 1 

is symmetric and positive definite, and /, g and ij: are determined by setting the exact solution 
u = sin(27r(t^ + 1) + 7r/2) sin(27rx + 7r/2) sin(27r?/ + 7r/2). The results are shown in Table 5 and 
Table 6, which indicate the optimal convergence rates for all norms presented. 



Table 5. Convergence rate for parabolic problem with inhomogeneous Dirichlet 
boundary condition with k = h 



h 


{oo,T} 


N^/i {oo.ar} 




{L2,r} 


\\^h\\{L'^^dT} 


1/8 


1.21E-01 


1.38E-01 


3.57E-01 


4.64E-02 


8.14E-02 


1/16 


5.64E-02 


6.10E-02 


1.23E-01 


1.87E-02 


3.39E-02 


1/32 


2.67E-02 


2.81E-02 


4.34E-02 


8.35E-03 


1.54E-02 


1/64 


1.29E-02 


1.33E-02 


1.55E-02 


3.95E-03 


7.29E-03 


1/128 


6.33E-03 


6.42E-03 


5.61E-03 


1.92E-03 


3.55E-03 


O(/i0r = 


1.0654 


1.1076 


1.4978 


1.1487 


1.1301 
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Table 6. Convergence rate for parabolic problem with inhomogeneous Dirichlet 
boundary condition with k = h? 



h 


\\^h\\{oo.T} 


\\^h\\{oo,dT} 


W^dehW 


{L2,r} 


{L2,aT} 


1/8 


7.41E-02 


9.68E-02 


1.24E-01 


3.53E-02 


5.74E-02 


1/16 


1.92E-02 


2.56E-02 


3.00E-02 


9.14E-03 


1.47E-02 


1/32 


4.83E-03 


6.44E-03 


7.45E-03 


2.30E-03 


3.70E-03 


1/64 


1.21E-03 


1.62E-03 


1.86E-03 


5.78E-04 


9.28E-04 


1/128 


3.35E-04 


4.34E-04 


4.67E-04 


1.53E-04 


2.48E-04 


0{K')r = 


1.9474 


1.9500 


2.0118 


1.9637 


1.9636 
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